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ABSTRACT 


The  problem  of  radiowave  propagation  over  irregular  terrain  is  solved  by  using  the 
standard  parabolic  equation  method.  The  ground  is  characterized  by  an  impedance  boundary 
condition  and  a  height  profile.  A  tropospheric  boundary  condition  is  used  to  truncate  the 
computational  domain.  This  thesis  uses  a  novel  approach  of  casting  all  the  equations  in  a 
curvilinear  coordinate  system.  The  coordinate  system  is  generated  in  a  simple  manner  using 
the  ground  profile  data.  The  equations  are  solved  by  the  finite  difference  method  using  the 
Crank-Nicolson  scheme. 

Different  numerical  values  for  various  important  parameters  (e.g.,  step  size,  location 
of  tropospheric  boundary,  the  region  above  the  tropospheric  boundary,  etc.)  were  used,  and 
their  effect  on  the  accuracy  and  computing  time  are  discussed.  Validation  of  the  numerical 
results  with  exact  and/or  experimental  results  are  presented  for  different  terrain  profiles.  Both 
perfectly  electric  conducting  (PEC)  and  lossy  impedance  surfaces  are  considered. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  topic  of  radiowave  propagation  over  irregular  terrain  is  extremely  important  in 
ground-to-ground  as  well  as  in  ground-to-air  communications  used  by  the  Navy,  Air  Force, 
and  the  Army.  Similarly,  the  ability  to  predict  radiowave  propagation  over  irregular  terrain 
has  a  significant  impact  in  determining  target  detectability  in  a  radar  system.  The  physics  of 
propagation  is  affected  by  ever-changing  atmospheric  conditions  and  by  complex  terrain 
features  on  the  ground.  The  link  reliability  in  a  communication  system  or  target  detectability 
in  radar  can  be  significantly  affected  by  the  so  called  'multipath  fading'.  The  path  between  a 
transmitter  and  a  receiver  is  often  obstructed  by  natural  or  man-made  obstacles  such  as  hills, 
buildings,  atmospheric  layers,  trees,  rain,  fog,  etc.  In  the  case  of  atmospheric  multipath  fading, 
abnormal  propagation  of  electromagnetic  waves  resulting  from  super-refraction  or  sub- 
reffaction  can  result  in  severe  loss  of  signals.  Reflection  multipath  fading,  which  is  due  to 
interference  between  the  direct  and  the  ground  reflected  waves  depends  strongly  on  the 
terrain  geometry  and  ground  constants.  It  is  therefore  important  to  designers  and  operators 
of  communications  and  radar  systems  to  predict  the  electromagnetic  fields  due  to  radiating 
sources  in  the  troposphere  and  to  assess  the  effects  of  environment  on  radiowave 
propagation. 

Any  in-depth  understanding  of  these  systems  requires  a  knowledge  of  the  physical 
phenomena  that  governs  low-angle  propagation;  more  specifically,  to  be  able  to  model 
complex  fading  phenomena  due  to  refiaction,  reflection,  scattering,  and  diffraction.  Numerous 
analytical  methods  are  available  for  predicting  electromagnetic  wave  propagation  such  as 
geometric  optics,  physical  optics,  normal  mode  analysis,  and  combinations  of  the  above. 
However,  they  have  several  limitations  for  predicting  propagation  in  a  complex  environment. 
The  parabolic  equation  method  has  been  used  to  predict  radiowave  propagation  in  an 
inhomogeneous  atmosphere  and  over  flat  terrain,  and  also  for  predicting  radiowave 
propagation  over  sloping  irregularities. 
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In  this  thesis,  we  will  concentrate  solely  on  the  effects  of  irregular  terrain.  Since  the 
propagation  path  could  extend  over  several  thousands  of  wavelengths,  it  is  important  to  have 
a  method  that  is  efficient  numerically.  The  parabolic  equation  method  is  one  such  method. 
One  advantage  of  a  parabolic  partial  differential  equation  (PDE)  over  an  elliptic  PDE  is  that 
in  the  former  case,  the  field  at  any  location  can  be  computed  in  terms  of  the  field  at  a  previous 
location.  However,  this  would  be  accurate  only  when  the  waves  propagate  predominantly  in 
the  forward  direction.  In  deriving  a  parabolic  PDE,  it  is  assumed  that  the  waves  are 
predominantly  forward  traveling.  This  is  approximately  met  in  a  typical  radio  link  where  the 
received  signal  is  primarily  affected  by  the  nature  of  the  path  between  the  transmitting  and 
receiving  antennas.  Of  course  there  are  several  situations  when  this  is  not  true.  For  example, 
when  there  is  a  large  obstacle  behind  a  receiving  anteima,  back-scattering  from  the  obstacle 
will  affect  the  received  signal.  Nonetheless,  the  parabolic  equation  method  has  been  used 
successfully  in  the  past  to  predict  propagation  in  several  scenarios,  particularly  for 
atmospheric  multipath  fading. 

B.  OBJECTIVE 

In  this  thesis,  we  adopt  the  same  parabolic  PDE  as  in  the  previous  approach  [Ref  1] 
to  predict  radiowave  propagation  over  an  irregular  terrain.  A  tropospheric  boundary  condition 
is  used  to  truncate  the  computational  domain  at  the  top,  while  an  impedance  boundary 
condition  is  used  on  the  uneven  terrain  to  characterize  the  ground.  The  key  feature  of  this 
thesis  is  to  use  a  novel  approach  of  casting  all  the  equations  in  a  curvilinear  coordinate 
system.  A  body  fitted  coordinate  system  is  generated  based  on  the  specification  of  the  ground 
profile.  This  permits  accurate  modeling  of  the  boundary  conditions  which  is  so  vital  to  the 
success  of  the  model.  The  parabolic  equation  method  is  a  full-wave  method  in  that  it  includes 
all  aspects  of  wave  propagation  such  as  forward  reflection,  refi-action,  diffraction,  and  surface 
wave  propagation.  However,  as  stated  above,  it  ignores  back-scattering.  Chapter  II  presents 
the  derivation  and  formulation  of  the  governing  partial  differential  equation  for  the  standard 
parabolic  equation  method,  the  impedance  boundary  condition,  and  the  tropospheric  boundary 
condition.  Also  presented  in  Chapter  11  is  the  transformation  of  all  the  partial  differential 
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equations  to  a  curvilinear  coordinate  system.  Chapter  III  details  the  generation  of  the 
coordinate  system  and  the  numerical  procedure  for  solving  the  parabolic  PDE.  The 
performance  of  the  numerical  solution  is  examined  in  Chapter  IV.  This  includes  a  study  on 
the  effects  of  using  different  numerical  values  for  various  important  parameters  (e  g.,  step 
size,  location  of  the  tropospheric  boundary,  the  region  above  the  tropospheric  boundary,  etc.) 
on  the  accuracy  of  the  solution,  and  validation  of  the  numerical  results  with  exact  and/or 
experimental  results  for  different  terrain  profiles.  Recommendations  and  conclusions  are 
presented  in  Chapter  V.  Finally,  MATLAB  computer  codes  for  the  numerical  implementation 
are  presented  in  the  Appendix. 
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n.  FORMULATION 


In  this  chapter  we  present  the  governing  partial  differential  equation  for  the  fields 
together  with  the  required  boundary  conditions.  We  present  only  the  final  forms  and  refer  the 
reader  to  [Ref  2]  for  more  details.  Figure  1  shows  a  Hertzian  electric  source  placed  over  an 
irregular,  lossy  terrain.  The  terrain  is  characterized  by  its  height  profile  and  an  impedance 
boundary  condition.  The  impedance  of  the  ground  depends  on  the  ground  constants  (e„ep 
PoPp  and  a).  We  wish  to  solve  the  fields  at  a  point  on/over  the  ground  in  the  presence  of  the 
irregularities.  We  consider  only  a  2-dimensional  situation  where  the  sources,  geometry,  and 
all  fields  are  z-invariant. 


Electric 


Source 


(xO,  yO) 


Figure  1 .  An  electric  source  producing  fields  over  an  irregular  terrain. 


Section  A  deals  with  the  impedance  boundary  condition.  In  Section  B,  we  present  the 
theory  on  the  parabolic  PDE.  In  Section  C,  we  present  the  tropospheric  boundary  condition 
required  to  terminate  the  computational  domain.  In  any  partial  differential  equation,  proper 
imposition  of  boundary  conditions  is  very  critical  to  the  final  solution.  We  desire  to  solve  the 
parabolic  PDE  by  an  implicit  finite  difference  scheme.  Most  of  the  previous  work  in  this  area 
used  a  cartesian  mesh  because  of  its  simplicity.  However,  practical  geometries  seldom 
conform  to  cartesian  coordinates.  Some  sort  of  interpolation  is  needed  near  the  boundary 
when  non-cartesian  geometries  are  encountered,  as  in  the  present  case.  This  could  result  in 
a  severe  loss  of  accuracy.  To  better  model  the  boundary  conditions,  we  solve  the  equations 
in  a  curvilinear  coordinate  system  generated  by  treating  the  lower  irregular  boundary  as  one 
coordinate  line.  In  Section  D,  we  cast  all  the  equations  in  a  curvilinear  coordinate  system. 

A.  IMPEDANCE  BOUNDARY  CONDITION 

An  impedance  boundary  condition  (IBC)  relates  the  tangential  components  of  the 
electric  and  magnetic  fields  at  the  interface  of  two  media.  If  v  is  a  unit  normal  vector  to  the 
boundary,  and  s  is  a  unit  tangential  vector  as  shown  in  Fig.  1,  the  boundary  condition  is  [Ref 
3] 

V  X  (  V  X  £  )  =  -ti^A  V  X  H,  (1) 

where  =  ZJr\^  is  the  surface  impedance  normalized  to  the  ffee-space  impedance 

actual  surface  impedance  that  is  dependent  on  the  media  constants 
and  the  incident  angles.  The  surface  impedance  is  determined  fi'om  the  intrinsic  impedance  of 
the  medium  by  considering  plane  wave  reflections  fi-om  the  interface.  Figure  2  shows  the 
interface  between  two  media.  The  complex  propagation  constants,  Yi,  Y2>  the  intrinsic 
impedances  rji  and  Tjj  can  be  expressed  in  terms  of  the  media  constants.  They  are 

Y?  (2) 

Y2  =  *  ju>ee)  =  (3) 
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Figure  2.  Boundary  interface  between  two  regions. 


and 


(4) 


According  to  Snell's  law,  Yi  sin  Oj  =  Y2  sin  0,.  The  plane-wave  reflection  coefficients  for  the 
vertical  and  horizontal  polarizations,  Ry  and  R^,  are  [Ref  4] 

„  iljsece,  -  n^sece, 
iljsece,  ♦  tioSecO, 


iljCosO,  -  11,COS0j 
lljCOS  0^  +  lljCOS  0, 


(6) 


Using  these,  the  surface  impedance  of  the  lossy  ground  can  be  taken  as 


Jf 

J  Horizontal  Polarization 

I  Z/'  =  T^^cos  0^,  Vertical  Polarization 


(7) 


Therefore, 


For  the  special  case  of  normal  incidence,  ij/j  =  90°, 


where  =  - .  In  this  thesis,  we  will  use  a  normalized  surface  impedance  given  in  Eq.  (9) 

for  all  the  resists.  For  a  2-Dimensional  case  with  vertical  polarization,  the  impedance 
boundary  condition  given  by  Eq.  (1)  can  be  simplified  as 

—  -  =0.  me  Vertical  Polarization  (10) 

Similarly,  we  have  for  the  horizontal  polarization 


1 

- Horizontal  Polarization 

dv 

“  jr 


(11) 


For  a  perfectly  conducting  material,  =  0  and  the  impedance  boundary  condition  reduces 
dH  ^ 

to  =  0  (Vertical  Polarization)  and  E^  =  Q  (Horizontal  Polarization). 

B.  STANDARD  PARABOLIC  PARTIAL  DIFFERENTIAL  EQUATION 

For  a  two-dimensional  electric  source  producing  fields  in  a  homogeneous  region,  all 
quantities  are  independent  of  the  z-coordinate,  and  propagation  takes  place  in  the  x^^-plane. 
From  Maxwell's  equations,  we  have  V  x  £  =  and  V  x  jy  =  yweE  >  J.  The  fields 

could  be  expressed  in  terms  of  the  z-component  of  the  magnetic  field,  in  the  case  of 
vertical  polarization  {TE^  fields),  and  in  terms  of  E^  for  the  horizontal  polarization  (TM^ 
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fields).  In  a  source-free  environment,  the  equation  satisfied  by  the  magnetic  field  is 

V  •  (VH^  *  =  0.  Vertical  Polarization  (12) 

Similarly  the  equation  for  the  electric  field  is 

V  •  (V£^)  ♦  klE^  -  0.  Horizontal  Polarization  (13) 

Equations  (12)  and  (13)  may  be  combined  into  an  equation  of  the  form 

V^ilr  .  =  0,  (14) 

where  for  TE  Polarization  and  ^  for  TM  Polarization.  Equation  (14)  is  the 

standard  Helmholtz  equation  and  is  elhptic  in  nature;  the  field  at  any  one  point  depends  on 
field  at  every  other  point  in  a  complicated  manner.  However,  for  wave  propagation  problems, 
an  approximate  answer  can  be  obtained  by  the  use  of  a  parabolic  PDE.  In  this  case,  the  field 
at  a  particular  range  depends  on  the  field  at  previous  range  points  only.  Assuming  that  the 
wave  propagates  predominantly  in  the  positive  x-direction,  we  write 

=  e*^u(xy),  (15) 


where  m(xj^)  is  a  slowly  varying  function  of  x.  We  now  impose  the  restriction  that 


l«J  «  I 


.  du  d^u 

dx  dx^ 


)  into  Eq.  (14)  and  arrive  at 

du  -J  d\ 
dx  2k ^  dy^ 


(16) 


(17) 


Equation  (17)  is  the  exact  form  of  the  narrow  angle  parabolic  PDE  approximation.  The 
impedance  boundary  conditions  derived  in  the  previous  section  can  also  be  expressed  in  terms 
of  the  V  function  : 

“v  ■  ° 


IBC  Vertical  Polarization 


(18) 


IBC  Horizontal  Polarization 


where 


—  =  -  sin  6  and  u 
dv 


du 

dv 


on  the  irregular  terrain.  By  defining 


-  sine) 
~  sine) 

I  Af 


Vertical  Polarization 
Horizontal  Polarization  ’ 


(19) 


(20) 


Eqs.  (18)  and  (19)  can  be  combined  and  written  as 

“v  *  =  0-  (21) 

The  parabolic  PDE  given  by  Eq.  (17)  is  valid  for  propagation  angles  close  to  the  horizontal 
(±15°  in  practice)  [Ref  5], 


C.  TROPOSPHERIC  BOUNDARY  CONDITION 

Our  computational  domain  consists  of  the  region  above  the  lossy  ground.  To  truncate 
the  computational  domain  on  the  upper  side,  we  consider  a  point  y  ~  yo  high  enough  and 
impose  a  boundary  condition  of  the  form 

"y  ♦  =  P  on  y  =  y^.  (22) 


To  derive  this,  we  start  with  the  parabolic  equation  -  2jku^  =  Owith  a  complex  k. 
Consider  the  problem  of  determining  the  field  at  any  point  in  x  >  x,„,  and  y>  y^,  given  the 
initial  data  on  x  =  x^„  m(x^,j;)  =J{y),  and  the  boundary  data  ony=y^  For 

analytical  simplicity,  we  assume  k  =  k^-]e,e>  0.  The  lossless  case  which  we  are  considering 
in  the  thesis  can  be  considered  as  the  limit  of  the  lossy  case  as  e  ->■  0.  Using  the  Fourier  sine 
integral,  it  is  shown  in  [Ref  2]  that 
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The  integrals  can  be  evaluated  approximately  by  replacing  the  derivatives  with 
differences.  Figure  3  shows  data  points  on  the  line  j  =yg  and  on  the  portion  of  x  =  x,„,  which 

is  above 


Figure  3.  Discrete  derivatives  for  integral  evaluations. 


Consider  the  evaluation  of  duldy\y=yo  at  x  -  x^„,  =  x^jj  with  initial  data  on  the  line  x  =  x,„^.  Let 
us  assume  that  this  initial  data  is  known  on  a  uniform  gridy„  =  Jo  +  wAj,  m  =  0,  1,  .  The 
derivative  in  the  interval  (y„.,,  j„)  can  be  approximated  by  the  forward  difference  formula 


aXy)  ^  -  “m-l 

dy  Ly 

where  u„  =  m(x,„  Then 


(24) 


—  f"" 


(y-- 


7l(x-xJ 


(y»-r 
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where  F(x)  =  is  a  complex  Fresnel  integral  [Ref.  6], 

Now,  consider  the  boundary  data  on  the  line  Assume  a  non-uniform  grid 

^ini  ^ini  ^2  •••»  ^ini  ^  On  the  interval  x  -  x„),  the  derivative  can  be 

approximated  as 


(26) 


where  —  u(x„,  y^.  The  second  integral  in  Eq.  (23)  can  now  be  obtained  as 


Substituting  Eq.  (25)  and  (27)  into  Eq.  (23), 


where  =  '/i(x^,+  xp  x^^-  x^,  =  Viix-  x^,)  &  'A  A  x^,.  Equation  (28)  can  be  thought 
of  as  the  discrete  version  of  a  continuous  boundary  condition  of  the  form 


where 


—  ♦  r(x)«  =  js:(x), 
dy 


(29) 
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D.  TRANSFORMATION  TO  A  CURVILINEAR  COORDINATE  SYSTEM 

The  partial  differential  equation  in  Eq.  (17)  and  the  boundary  conditions  in  Eqs.  (21) 
and  (28)  will  now  be  transformed  to  a  curvilinear  coordinate  system.  Consider  the  narrow 
angle  parabolic  PDE  given  by  Eq.  (17) 


u 


X 


1  d^u 


(32) 


together  with  the  tropospheric  boundary  condition 

— (33) 
dy 

and  the  impedance  boundary  condition  on  the  irregular  boundary 

♦  Cjtt  =  0.  (34) 

We  will  cast  all  of  these  equations  in  a  curvilinear  coordinate  system  (^,  q)  shown  in  Fig.  4. 
This  will  permit  accurate  imposition  of  boundary  conditions.  Note  that  the  parabolic  nature 
of  the  equation  will  not  change  as  a  result  of  the  transformation.  The  irregular  terrain 
boundary  will  map  into  tj  =  0  and  the  upper  boundary  into  N  (integer).  The  entire 
coordinate  system  is  generated  by  the  specification  of  the  terrain  geometry. 

Assuming  the  transformation  x  =  x(Q,  y  =  K^.ti),  and  using  the  various  metrics 
needed  [Ref  7],  Eq.  (32)  becomes 
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Letting 


Boundaiy  Condition  2 

Physical  Domain  Computational  Domain 

Figure  4.  Curvilinear  coordinate  system. 


then  Eq.  (35)  can  be  expressed  as 


b^U 


’ll  ■ 


The  normal  derivative  on  a  t)  =  constant  line  is 


u  (Tf^constant  )  = 


2  2 


~u  - 


I 

Vi  >'{*1  ^xl  .  yl  (»jV^  - 


(39) 


Using 


cos  6  = 


and  sin  6 


(40) 


f~2  2  ri  2 

Vi  *  yi  vi  * 

the  boundary  condition  at  the  bottom  boundary  +  C;«  =  0  gets  transformed  for  =  0  to 


Tl 

-  — sinOcosBw,  ♦  cos  6u  =  0. 

Similarly,  the  boundary  condition  at  j  gets  transformed  to 

«,(£..  S)  .  r(5.,  JV)y,(£>(E„  J^)  .  y,J(£.,  AT). 


(41) 


(42) 
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m.  SOLUTION  PROCEDURE 


In  the  previous  chapter  we  presented  the  underlying  equations  and  the  relevant 
boundary  conditions  for  the  computation  of  fields  over  an  irregular  terrain.  In  Section  A  of 
this  chapter,  we  present  details  on  the  generation  of  the  coordinate  system.  In  Section  B,  we 
adopt  the  Crank-Nicolson  scheme  in  the  transformed  domain  to  obtain  the  required  finite 
difference  equations. 

A.  GENERATION  OF  THE  CURVILINEAR  COORDINATE  SYSTEM 

Consider  the  computational  domain  between  the  terrain  profile  and  a  horizontal  upper 
boundary.  We  treat  the  terrain  profile  as  being  comprised  of  piece-wise  linear  functions  and 
consider  a  horizontal  upper  boundary  as  shown  in  Fig.  5.  We  wish  to  generate  a  coordinate 
system  such  that  the  lower  and  the  upper  boundaries  correspond  to  constant  r\  lines. 
The  coordinate  lines  gradually  become  horizontal  as  one  proceeds  from  the  lower  to  the 
upper  boundary.  Spedfication  of  the  lower  and  upper  boundaries  completely  determines  our 
coordinate  system.  Without  loss  of  generality,  we  will  choose  the  increments  between  the 
mesh  points  in  the  (^,11)  plane  to  be  unity.  Figure  5(a)  shows  the  mesh  points  in  the  physical 
domain  and  Fig.  5(b)  in  the  computational  dommn.  Figure  6  shows  corresponding  pairs  of 
points  on  the  upper  and  lower  boundaries  in  the  physical  and  computational  domains.  The  (x, 
y)  coordinates  of  an  interior  point  are  generated  by  using  linear  interpolation.  Letting  Ax,  = 
(x,+/  -  x^.  Ay,  =  (v,+/  -y^,  and  A  5i  =  (^i+i  -  Q,  the  parametric  equations  for  x  andy  are  [Ref 
2] 

xil)  -  X,  .  ^(5  -  1),  (43) 

yil.Tt)  .  1 1  -  ij[  y,  .  ^(5  -  g]  .  iy..  (44) 

for  0  <  Ti  i  N. 
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.(a) 

ilcal 


Physical  Domain 


(b) 

Computational  Domain 


Figure  5.  Physical  and  computational  grids. 


Figure  6.  Physical  and  computational  domain  segments. 


It  is  possible  to  generalize  this  to  any  ground  profile  shape  of  the  form  f(x)  by  writing 

with  xiX)  given  as  above.  However,  we  choose  the  piece- wise  linear  form  for  its  ability  to 
model  discrete  sets  of  data  on  the  ground.  At  any  interior  point,  <  ^  <  the  various 
metrics  are  obtained  as 


At  the  junction  points  ^  or  5,+/, 


we  use  the  central  difference  formulas 


-  ** 

Ax,. 

1  " 

5..  -  £.• 

< 

1 

+  Ay, 

V'n. 

’a  5,. 

*  A5,  • 

(48) 


(49) 


Note  that  the  analytical  expression  yields  a  discontinuous  value  for  these  derivatives. 
We  use  the  analytical  expressions  only  to  generate  grid  points  and  use  the  central  difference 
formulas  to  arrive  at  the  derivatives  with  respect  to  In  other  words,  once  the  grid  points 
are  generated  on  the  lines  BR, ... ,  etc.,  we  assume  that  the  space  is  smoothly  connected 
through  the  grid  points.  In  the  numerical  implementation  using  the  Crank-Nicolson  implicit 
scheme  [Ref.  8],  the  metrics  are  needed  at  the  midpoint  with  respect  to  i.e.,  at 
A  5/2,  and  the  interior  point  formulas  are  applicable.  For  a  uniform  mesh  in  the  computational 
domain,  A  =  1,  t)  =q,q=^0,  1,2, ... ,  N. 
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=  Ax,, 

(50) 

< 

1 

II 

(51) 

> 

(52) 

q 

(53) 

2  )' 

—y  . 

N  ^ 

Figure  7  shows  an  example  of  the  curvilinear  mesh  for  a  Gaussian  shaped  ridge  on  flat 
ground.  Note  that  the  vertical  increment  is  constant  on  any  vertical  line,  but  varies  from  line 
to  line. 


Figure  7.  Curvilinear  mesh  for  a  Gaussian  shaped  ridge  on  flat  ground. 
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B.  NUMERICAL  IMPLEMENTATION  USING  THE  CRANK-NICOLSON 
SCHEME 

Consider  the  narrow  angle  parabolic  PDE  together  with  the  boundary  conditions  given 
by  Eqs.  (38),  (41),  and  (42).  We  would  like  to  implement  the  equations  using  a  Crank- 
Mcolson  implidt  scheme  [Ref.  8].  We  will  write  different  forms  for  the  internal  and  boundary 
regions.  The  grid  points  for  an  internal  region  in  the  computational  domain  are  shown  in  Fig. 
8.  Given  the  field  data  on  the  line p-\,  we  would  like  to  predict  the  field  on  the  line p. 


Figure  8.  Internal  grid  points  in  computational  domain. 

Using  the  notation  =  u{p,q)  =  u{x^^,  the  various  derivatives  (assuming  =  Aq  =  1) 
are 

«{(p-*/2,?)  =  «/  -  (54) 
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■  y.[<,-  2.,/.  2«;‘.  <■]. 


Substituting  into  Eq.  (43)  and  rearranging  the  terms,  we  get  [Ref.  2] 


Letting 


P  =  (*3)7^ 


Y  = 


(55) 

(56) 


(57) 

(58) 

(59) 

(60) 


we  have  for p=\,2,  •••;  ^=1,2,  A^-1 

“«/i  -  (1  -  P)«/  ♦  Y«/i  =  -  (1  -  P)«f'  -  y«7i'.  (61) 

We  will  extend  the  applicability  of  this  equation  to  include  q  =  0  and  q  =  N\.o  accommodate 
the  derivative  boundary  conditions  on  the  lower  and  upper  boundaries.  For  ^  =  0,  we  have 

au^  -  (1  +  p)u/  +  Y“.f  =  -aaf'  -  (1  -  P)^^'  -  Y«.7*-  (62) 

For  q  =  N,vfe  have 

-  (1  "  P)<  *  Y<i  =  -a<;  -  (1  -  P)«r  -  Y<;  (63) 
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Figure  9.  Grid  points  at  the  lower  boundary. 

Figure  9  shows  the  grid  points  at  the  lower  boundary.  We  use  central  differences  for  the 
boundary  condition  in  Eq.  (41)  and  combine  with  Eq.  (62)  to  get 

a'tt'  -  (1  .  -  (1  -  (64) 

where 

a'  =  a  +  Y,  (65) 

P'  =  P  -  ’  (66) 

p"  =  p  -  2Yy^cos0^c,  ♦ 

We  have  also  tried  to  use  one-sided  forward  difference  formulas  near  the  terrain,  but  the 
results  were  less  accurate  than  with  the  central  difference  scheme.  Results  for  the  two 
approaches  will  be  compared  in  the  next  chapter. 
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Figure  10.  Grid  points  at  the  upper  boundaiy. 


Similarly,  using  central  differences  for  the  tropospheric  boundary  condition  in  Eq.  (42)  and 
combining  with  Eq.  (63)  results  in 


(1  .  - 


=  (1  -  >«)« 


N 


yu 


p-i 

N-l  * 


(68) 


where 


A,  =  6  +  y  8a 


Jk 


71  Ax 


>1 


y  =  y  *  (c. 


'  ^7-  ■  Vi’  Vk  =  Vi  *  *  Vi>’ 
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We  augment  Eq.  (61)  for  q=  \,2,  ••• ,  N-\  with  Eq.  (64)  and  (68)  for  ^  =  0  and  q  =  N 
respectively  to  define  the  equations  for  ^  =  0,  1,  2,  N.  The  system  of  equations  so  defined 
can  be  expressed  as  a  matrix  equation  of  the  form 


X  X 
XXX 
XXX 


X  X 
XXX 
XXX 


4ay 


(69) 


where  ^  denotes  a  non-zero  entry.  The  tri-diagonal  matrix  on  the  left  hand  side  ofEq.  (69) 
can  be  inverted  efficiently  to  yield  a  solution  on  line  ^  in  terms  of  the  field  values  on  the 
line  ^  Equation  (69)  can  be  used  to  march  forward  in  range  starting  fi'om  initial  data 
specified  at  ^  =  0  (x  =  x,„,). 
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IV.  COMPARISON  OF  RESULTS 


Ha\dng  presented  the  basic  theory  in  the  previous  chapters,  we  will  present  numerical 
results  in  this  chapter.  The  overall  accuracy  of  the  solution  depends  on  a  number  of 
parameters  such  as  the  step  size,  location  of  the  tropospheric  boundary  (i.e.,  in  Fig.  1),  the 
region  above  the  tropospheric  boundary,  and  the  differencing  scheme  used  at  the  boundaries. 
In  this  chapter  we  will  make  a  systematic  study  on  the  dependence  of  the  solution  on  various 
parameters.  Section  A.  1  compares  the  results  of  the  forward  versus  central  differencing 
schemes  for  implementing  the  lower  boundary  condition.  Sections  A.2  and  A.3  show  the 
effects  of  horizontal  and  vertical  step  sizes,  respectively.  In  Section  A.4,  we  present  the  effect 
of  the  placement  of  the  tropospheric  boundary  condition.  In  the  remaining  sections,  we 
present  comparisons  with  the  results  avmlable  in  the  literature  for  specific  obstacles. 

A.  PARAMETRIC  STUDY  OF  THE  FINITE  DIFFERENCE  SCHEME 

1.  Forward-Difference  versus  Central  Difference  Scheme  for  the  Impedance 
Boundary  Condition 

In  the  implementation  using  the  Crank-Nicolson  implicit  scheme,  the  derivatives  at 
any  interior  point  are  approximated  using  the  central  difference  formulas.  This  approximation 
is  also  extended  to  the  lower  boundary  given  in  Eq.  (69).  Alternatively,  the  derivatives  at  the 
lower  boundaiy  could  also  be  approximated  using  the  forward  difference  formulas.  From  Eq. 
(43)  and  using  the  forward  difference  formulas,  we  have 


Equation  (70)  can  be  written  as 

-  (1  .  *  (1  -  PX""'.  (71) 
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where 


a 


1, 


ft,  J  2sin0 

P  =  y^cos  01  - 

\  Xr 


-4 


o  //  qi  2siii  0 

P  =  y^cos  0^ -  +  Cj 


(72) 


Numerical  results  using  both  the  central  difference  and  forward  difference  schemes 
were  used  to  compute  magnetic  field  due  to  a  vertically  polarized  line  source  on  the  surface 
of  a  PEC  plane.  The  results  of  both  schemes  are  compared  with  the  exact  solution  [Ref  9] 
and  plotted  in  Fig.  1 1  (magnitude  of  magnetic  field,  Hz,  on  the  surface  of  the  PEC  plane), 
Fig.  12  {Hz  @  Height  =  0.  U),  and  Fig.  13  {Hz  @  Height  =  20A).  It  is  observed  that  the 
central  difference  scheme  produces  a  solution  that  agrees  better  with  the  exact  solution.  The 
solution  using  forward  differences  exhibits  spurious  oscillations  which  could  not  be  resolved. 
Henceforth  we  will  adopt  the  central  difference  scheme  for  the  subsequent  test  cases  shown. 

2.  Variation  in  Horizontal  Step  Size  (With  Fixed  Vertical  Step  Size) 

The  Crank-Nicolson  scheme  has  the  advantage  of  being  valid  (i.e.,  convergent  and 
stable)  for  all  finite  values  of  the  horizontal  step  size  for  a  cartesian  mesh.  To  verify  the 
validity  of  the  scheme  with  respect  to  horizontal  step  sizes,  numerical  results  were  evaluated 
for  three  different  horizontal  step  sizes  (Ax  =  0.  U,  Ax  =  0.2  A,  and  Ax  =  0.5  A)  with  a  fixed 
vertical  step  size  of  Ay  =  0. 1  A.  The  magnitude  of  the  magnetic  field,  Hz,  due  to  a  vertically 
polarized  line  source  placed  at  {xO  =  0,y0  =  0)  on  the  surface  of  a  PEC  plane  are  computed 
numerically  and  plotted  in  Fig.  14  {Hz  on  the  surface  of  the  PEC  plane)  and  Fig.  15  {Hz  @ 
Height  =y„  =  20A).  The  plots  show  that  the  results  are  convergent  for  all  the  three  horizontal 
step  sizes,  and  excellent  agreement  is  observed  between  the  numerical  results  and  the  exact 
solution  for  the  surface  magnetic  field.  Minor  oscillations  are  observed  for  magnitude  of  the 
magnetic  field  at  the  upper  boundary.  However,  the  average  value  of  the  magnetic  field,  Hz, 
is  still  close  to  the  exact  solution  given  by  [Ref  9] 


H,{xy) 


jK  -  V 

4  R 

O 


(73) 
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where  x  =  kx  ,  =  kx^  ,  -  x^^  +  (y  -  y^,  and  is  the  Hankel 

function  of  the  second  kind  of  order  1. 

3.  Variation  in  Vertical  Step  Size  (With  Fixed  Horizontal  Step  Size) 

Next,  we  investigate  the  behavior  of  the  Crank-Nicolson  scheme  with  respect  to 

vertical  step  size.  In  this  study,  the  numerical  results  were  evaluated  for  three  different  vertical 
step  sizes  {Ly  =  0. 1  A,  A>>  =  0.2A,,  and  =  0.5X)  with  a  fixed  horizontal  step  size  of  Ax  = 
0.1  A,.  The  same  profile  of  a  vertically  polarized  line  source  placed  at  (xO  =  0,y0  =  0)  on  the 
surface  of  a  PEC  plane  is  adopted,  and  the  magnitude  of  the  magnetic  field,  Hz,  are  computed 
numerically.  The  results  are  plotted  in  Fig.  16  {Hz  on  the  surface  of  the  PEC  plane)  and  Fig. 
17  {Hz  @  Height  =  yo  —  20 A,).  Again  the  plots  show  that  the  results  at  the  surface  are 
convergent  for  all  the  three  vertical  step  sizes  tested,  and  excellent  agreement  is  observed 
between  the  numerical  results  and  the  exact  solution  for  the  magnetic  field,  Hz. 

Minor  oscillations  are  still  observed  at  the  upper  boundary.  Poor  result's  are  obtained 
with  A^  =  0.5X.  This  is  possibly  due  to  the  inaccuracies  of  the  upper  boundary  condition. 
Nonetheless,  the  numerical  computation  using  step  size  of  Ax  =  A>'  =  0. 1  A,  produces  results 
which  are  in  good  agreement  with  the  exact  solution. 

4.  Variation  in  the  Tropospheric  Boundary  Condition  Parameter 

In  the  evaluation  of  duldy\y=yg  at  x  =  x^,4+x,„„  we  have  considered  and  assumed  that 
initial  data  is  known  on  a  uniform  grid>'„  -  Jo  “  w  =  0, 1,  •••.  In  practice  we  stop  at  some 
height  >>^o.  To  ascertain  how  high  a  point  should  be  considered,  we  examine  the  cases 
of  extending  beyond  the  tropospheric  boundary  by  5A,,  lOA,,  and  20X,  (\.t.,y,j^=yo  +  5A,, 
10  A.,  and  20  A.).  The  simple  test  profile  of  the  propagation  of  the  magnetic  field,  Hz,  due  to 
a  vertically  polarized  line  source  placed  at  {xO  =  0,y0  =  0)  on  the  surface  of  a  PEC  plane  are 
computed  numerically.  In  this  study  we  choose  Ax  =  A^'  =  0. 1  A.  and  the  upper  boundary,  y^ 
=  20X.  The  results  are  plotted  in  Fig.  18  {Hz  on  the  surface  of  the  PEC  plane)  and  Fig.  19  {Hz 
@  Height  =  20X).  The  plots  show  that  the  results  are  almost  identical  for  all  the  three 
cases.  To  get  an  accurate  representation  of  the  tropospheric  boundary  condition,  we  need 
only  to  consider  data  on  a  vertical  line  which  is  within  5  A,  from  the  upper  boundary.  As  such, 
all  subsequent  test  cases  vrill  be  computed  using  a  tropospheric  boundary  with 
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Figure  1 2.  Hz  versus  distance  @  height  =  0. 1  A. 
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Figure  15.  Hz  versus  distance  @  height  =  y„  =  20 X. 
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Figure  1 8.  Surface  magnetic  field  versus  distance 


B.  PROPAGATION  OVER  A  PEC  PLANE 


In  this  section,  we  examine  the  effect  of  the  height  of  the  upper  boundary,  on  the 
stability  and  accuracy  of  the  numerical  results.  The  magnitude  of  the  magnetic  field,  Hz,  due 
to  a  vertically  polarized  line  source  placed  at  {xO  =  0,y0  =  0)  on  the  surface  of  a  PEC  plane 
are  computed  numerically.  In  this  study  we  choose  Ax  =  Aj  =  0. 1 A  and  y,^^  =y^  +  51. 
Results  are  presented  for  5  A,  10 A,  20 A,  and  30 A.  The  initial  data  line  is  at  =  5A,  and 
the  magnetic  field  is  marched  over  a  horizontal  distance  of  10 A. 

1.  Magnetic  Field  Variation  with  Horizontal  Distance 

The  variation  in  the  magnitude  of  the  magnetic  field,  Hz,  along  the  horizontal  distance 
on  the  surface  of  the  PEC  plane  and  at  a  height  of  5  A  is  plotted  in  Fig.  20  and  Fig.  21 
respectively.  We  noticed  that  the  oscillations  are  predominant  when^’^  is  less  than  20 A.  For 
y„  >  20  A,  the  results  are  in  good  agreement  with  the  exact  solution. 

2.  Magnetic  Field  Variation  with  Vertical  Distance 

Figures  22,  23,  and  24  depict  the  results  of  the  variations  of  the  magnitude  of  the 
magnetic  field,  Hz,  with  vertical  height  at  horizontal  distances  of  6A,  lOA,  and  15A  from  the 
transmitter.  Again,  oscillations  are  present  when  the  height  of  the  upper  boundary  is  below 
20 A.  More  importantly,  we  observed  that  for  a  given  horizontal  distance,  the  deviations  of 
the  numerical  results  from  the  exact  solution  become  more  pronounced  at  higher  points.  This 
phenomenon  is  expected  because  the  parabolic  equation  used  in  this  thesis  is  valid  for 
propagation  angles  close  to  horizontal  (±15°).  In  fact  the  numerical  result  is  fairly  consistent 
with  this  restriction.  For  example,  at  a  horizontal  distance  of  15  A  fi'om  the  transmitter  (i.e., 
lOA  from  the  initial  data  line),  the  maximum  height  for  which  good  results  are  expected  is 
between  0  and  1.763  A.  From  Fig.  24,  we  see  that  this  is  true,  and  deviation  between  the 
numerical  results  and  the  exact  solution  begins  at  approximately  1.7A.  The  same  trend  is  also 
observed  in  Fig.  22  and  23.  However,  the  difference  between  the  numerical  results  and  the 
exact  solution  increases  with  distance  fi'om  the  transmitter.  This  could  be  due  to  inaccuracies 
of  the  tropospheric  boundary  condition.  The  error  in  the  tropospheric  boundary  condition  is 
accumulative,  and  we  would  expect  to  see  higher  error  at  larger  distances. 
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Figure  23.  Hz  versus  height  @  distance  =  10 A,. 


Figure  24.  Hz  versus  height  @  distance  =  1 5  A. 


C.  PROPAGATION  OVER  A  PEC  KNIFE-EDGE 

Consider  a  perfectly  conducting  knife  edge  of  height  as  shown  in  Fig.  25,  set  on  a 
perfectly  conducting  plane,  with  the  transmitter  and  receiver  both  on  the  ground  at  distances 
d,  and  The  attenuation  factor,  AF,  at  point  B  can  be  found  analytically  and  is  given  by 
Plef  10] 


AF  =  4 


<4i, 

2  Ju 


where 


«(0  =  i 


2d 


“c  -  “(5=*c)  =  K 
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2d 


(74) 


(75) 


Figure  25.  Perfectly  conducting  knife  edge  between  the  transmitter  at  A  and  the  receiver 
at  B,  both  of  which  are  on  a  perfectly  conducting  ground. 
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We  consider  the  case  where  =  2X  and  d,  =  \0X.  For  this  choice  Sj  =  1 1.3°.  A  vertically 
polarized  line  source  is  placed  at  x0^0,y0  =  0  and  the  magnetic  field  is  determined  from  an 
initial  range  of  5  A  to  a  distance  of30A.  The  height  of  the  upper  boundary  is  =  20A  andy^f,^ 
=To  +  5A. 

The  variation  of  the  surface  magnetic  field  over  the  horizontal  distance  is  plotted  in 
Fig.  26.  Variation  of  the  magnetic  field  with  height  at  a  distance  =  2k  is  plotted  in  Fig.  27. 
Comparison  of  the  attenuation  factors  at  d2  =  2k,  lOA,  15A,  and  20k  (02  =  45°,  11°,  7.6°, 
5.7°  respectively)  computed  numerically  and  the  exact  solution  using  Eq.  (74)  is  given  in 
Table  1. 


2k 

0.7758 

0.4071 

90.57% 

lOA 

1.1674 

0.6725 

73.59% 

15A 

0.9330 

0.7250 

28.69% 

20A 

0.8048 

0.7560 

6.46% 

Table  1 .  Attenuation  factor  for  propagation  over  PEC  knife  edge. 

The  rather  large  oscillations  seen  near  the  surface  in  Fig.  27  could  be  due  to  the  non¬ 
applicability  of  the  parabolic  equation  for  high  angles  (i.e.,  02  >  15°). 

D.  PROPAGATION  OVER  A  CIRCULAR  BOSS  IN  A  PEC  PLANE 

Next,  we  compare  the  numerical  results  for  propagation  over  a  semi-circular  boss  on 
a  perfectly  conducting  plane.  This  problem  typifies  radiowave  propagation  over  uneven 
terrain.  A  vertically  polarized  line  source  is  placed  in  front  of  the  semi-circular  boss  of  radius 
b  =  0.5  k  and  centered  at  the  origin.  The  transmitter  is  at  xO  =  -0.75  A,  y0  =  0.lk  from  the 
origin  and  the  initial  data  line  is  at  -0.6A  from  the  center  of  the  boss.  There  is  good  agreement 
between  the  exact  solution  given  in  [Ref  9]  and  the  numerical  results  generated  here  as 
evidenced  in  Fig.  28. 
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Figure  26.  Normalized  surface  magnetic  field  versus  distance  (PEC  knife  edge). 
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Figure  27.  Normalized  magnetic  field  versus  height  @  d2  =  2X  (PEC  knife  edge). 
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Figure  28.  Magnitude  of  the  normalized  surface  fields  versus  horizontal  distance  for 
a  line  source  placed  over  a  perfectly  conducting  plane  with  a  semi-circular  boss. 
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E.  PROPAGATION  OVER  A  LOSSY  IMPEDANCE  PLANE 

In  the  previous  sections,  we  have  considered  only  propagation  over  PEC  surfaces.  In 
this  section,  we  will  investigate  the  case  of  propagation  over  a  lossy  impedance  plane.  The 
ground  constant  chosen  are  e,  =  10,  o,  =  180  (corresponding  to  a  =  10  mS/m  at  1  MHz). 
The  line  source  is  placed  at  xO  =  0,y0  =  0.01  A.  The  horizontal  step  size.  Ax  =  0. 1  A,  and  the 
vertical  step  size,  Aj  =  0. 1  A.  The  initial  data  line  is  5  A  away  from  the  source,  and  the  solution 
is  allowed  to  propagate  for  25  wavelengths.  The  magnetic  field  on  the  initial  data  line  i 
generated  from  the  fields  on  a  flat  lossy  plane  [Ref  11] 


is 


H  = 

J 


where 


k 
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COS0^®(A:/,)  ♦  -  J2A^ 
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.  vV: 
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"i  =  \/(x  -  xf  *{y  -  yf  ,  ^(x  -  xf  .  (y  *  yf  , 


cos  6, 


cos  02  = 


Oc  - 


sin  02  = 
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The  results  for  the  surface  magnetic  field  are  plotted  in  Fig.  29.  Comparison  is  made 
with  the  data  points  for  the  exact  solution  given  in  Fig.  22,  pp.  52  of  [Ref  9].  There  is  very 
good  agreement  between  the  numerical  results  and  the  exact  solution. 


F.  PROPAGATION  OVER  A  LOSSY  GAUSSIAN  Hn  .T. 

Next,  we  treat  the  case  of  propagation  over  a  Gaussian  shaped  ridge  or  hill.  The 
height  h(x)  of  the  Gaussian  hill  is  defined  as 

h(x)  =  He  (77) 

where  H  is  the  maximum  height,  c  is  the  location  of  the  maximum,  and  w  is  a  parameter 
controlling  the  width  of  the  hill.  A  vertically  polarized  line  source  was  placed  at  xO  ^0,y0  = 
0.01  A  with  H  =  \  km,  c  =  5  km,  and  w  =  3  km.  The  initial  data  line  was  at  x,„,  =  2  km. 
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Figure  29.  Surface  magnetic  field  versus  distance  (lossy  plane). 
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Figure  30.  Magnitude  of  the  normalized  surface  fields  versus  horizontal  distance  for 
source  placed  at  the  origin  in  front  of  a  Gaussian  hill. 


The  ground  constants  were  6,=  10  and  a,  =  180  (corresponding  to  a  =  10  mS/m  at  1  MHz). 
We  plot  the  results  of  the  normalized  magnetic  field  in  Fig.  30  together  with  the  results 
generated  from  WAGNER  [Ref  12].  Although  there  is  minor  differences  in  the  magnitude 
of  the  normalized  magnetic  field,  it  is  seen  that  the  trend  in  the  two  results  are  in  good 
agreement.  Both  results  predict  large  peak  around  3.9  km  fi’om  the  source,  which  is  attributed 
to  focusing  by  the  concave  surface  of  the  hill.  A  secondary  peak  at  4.5  km  is  also  predicted 
by  the  numerical  method. 

G.  PROPAGATION  OVER  A  PEC  ISOSCELES  HILL 

Experimental  results  for  propagation  over  a  PEC  isosceles  triangular  hill  are  given  in 
[Ref  10].  In  [Ref  10],  an  irregular  path  consisting  of  an  isosceles  hill  of  height  1.23X  was 
constructed  fi-om  aluminum  sheets  whose  surface  impedance  was  taken  to  be  zero.  The 
dimensions  of  the  sloping  side  of  the  hill  are  shown  in  Fig.  31. 


Initial  Data  Line  Hill  H  =  3.85  cm 


Figure  3 1 .  Perfectly  conducting  isosceles  triangular  hill  of  height  1 .23  X  and  baselength 
20.16A,. 
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Figure  32.  Magnitude  of  the  normalized  fields  across  a  perfectly  conducting  isosceles 
triangular  hill. 
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Using  the  vertical  line  sources  at  a  distance  of  10.24  A  from  the  beginning  of  the  hill, 
and  a  frequency  of  9.6  GHz,  the  numerical  solution  is  computed  and  marched  from  an  initial 
data  line  of  5. 12 A.  The  normalized  magnetic  field  taken  at  a  height  =  0.16A  above  the 
irregular  surface  is  plotted  in  Fig.  32.  We  have  computed  the  attenuation  factor  for  the 
experimental  values  given  in  Fig.  4,  pp.  39  of  [Ref  10]  as  follows 


AFidB)  =  20  log 


10  log 


2ind  compared  them  to  the  results  obtained  from  our  numerical  method.  In  our  case,  we  have 
chosen  r,  =  5. 12A  to  coincide  with  the  initial  data  line.  It  can  be  seen  that  there  is  excellent 
agreement  between  the  experimental  data  and  the  results  from  the  numerical  solution. 

H.  PROPAGATION  OVER  A  PEC  CLIFF 

Measured  data  for  the  magnitude  of  the  field  strength  over  an  aluminum  cliff  edge  was 
also  given  in  [Ref  10].  Measurements  were  carried  out  over  the  aluminum  cliff  which  was 
situated  10.24A  from  the  transmitter  and  of  height  1.05A  as  shown  in  Fig.  33. 


Transmitter  Initial  Data  Line 


uiq>erturbed  groimd 


3.27  cm 
1.05X) 
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We  set  up  our  computation  model  based  on  the  geometry  given  in  Fig.  33.  At  a 
frequency  =  9.6  GHz,  we  placed  a  vertical  line  source  at  10.24A  from  the  edge  of  the  cliff, 
and  set  the  initial  data  line  at  5.12X.  We  have  chosen  Ax  =  Ay  =  0.08  A  (instead  of  Ax  =  Ay 
—  0. 1  A,)  to  enable  plotting  of  data  at  a  height,  =  0. 16A  (2  vertical  steps)  at  which  height  the 
measurements  were  also  taken.  The  result  of  the  numerical  solution  is  plotted  in  Fig.  34.  From 
the  plot,  we  observe  that  the  field  attenuates  over  the  first  flat  section  as  if  the  cliff  is  were 
absent.  Immediately  beyond  the  cliff,  however,  the  field  strength  drops  suddenly  in  the 
shadow  region,  and  then  increases  with  further  increase  in  range,  before  returning  as  before 
to  decreasing  values  with  increased  range.  We  have  computed  the  attenuation  factor  based 
on  the  measured  data  given  in  Fig.  6,  pp.  39  of  [Ref  10]  using  Eq.  (78)  and  superimposed 
them  in  Fig.  34.  Again,  there  is  excellent  agreement  between  the  numerical  solution  and  the 
experimental  results,  even  for  points  near  the  cliff  edge  and  immediately  beyond  the  cliff. 

I.  PROPAGATION  OVER  A  DIELECTRIC  COATED  CLIFF 

Experimental  results  were  also  available  for  a  dielectric  coated  cliff  [Ref  10]. 
Measurement  have  been  carried  out  over  the  same  cliff  edge  given  in  Fig.  33,  but  with  the 
horizontal  surfaces  constructed  from  polypropylene  covered  aluminum  sheet  and  the  vertical 
edge  made  of  aluminum  sheet  alone.  The  normalized  surface  impedance  of  the  polypropylene- 
covered  sheet  with  e,  =  2.2  and  dielectric  thickness,  can  be  calculated  to  be  A,”  =  jO.  171, 
which  is  highly  inductive.  We  use  the  same  computational  model  as  in  Section  H  and  plot  the 
results  in  Fig.  35.  Again,  we  compute  the  attenuation  factor  based  on  the  measured  data  given 
in  Fig.  7,  pp.  40  of  [Ref  10]  using  Eq.  (78)  and  superimposed  it  in  Fig.  35.  Again,  we  observe 
very  good  agreement  between  the  magnitude  of  the  of  the  attenuation  factor  between  the 
numerical  solution  and  the  measured  value.  This  is  true  for  points  near  the  cliff  edge  and 
immediately  beyond  the  cliff  edge.  The  behavior  of  the  field  strength  computed  numerically 
also  corresponds  very  well  with  the  experimental  results.  The  results  in  Sections  H  and  I 
clearly  demonstrated  the  versatility  of  this  numerical  method  for  terrain  with  PEC  surface  and 
dielectric  coated  surface. 


56 


Attenuation  Factor  (dB) 


Figure  34.  Magnitude  of  the  normalized  fields  across  a  perfectly  conducting  cliff. 
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Figure  35.  Magnitude  of  the  normalized  fields  across  dielectric  coated  horizontal 
surfaces  on  a  cliff  edge;  =  2.2.  Dielectric  thickness  0.15  cm.  Normalized  surface 
impedance  is  jO.  1 7 1 . 
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J.  PROPAGATION  OVER  A  MIXED  PATH  (LAND-SEA-LAND) 

Next,  we  would  like  to  examine  the  capability  of  the  parabolic  equation  formulation 
to  predict  radiowave  propagation  over  long  range  and  inhomogeneous  ground.  Consider  a 
four-thirds  earth  atmosphere  over  a  ground  with  electrical  properties  that  change  along  the 
propagation  path.  The  initial  and  final  portions  of  this  "mixed-path"  problem  are  over  land  (e, 
=  15,  a  =  0.03  mho/m),  while  fi’om  47  to  50  km  fi’om  the  source  the  propagation  takes  place 
over  sea  (e,  =  80,  a  =  4  mho/m).  Results  of  path-loss  predictions  at  a  receiver  height  z,.  =  5 
m  for  a  59.7  MHz  transmitter  at  a  height  of  z,  =  5  m  using  BFDG  method  with>’„  =  10m, 

=  40  km,  =  Im,  and  Ax  =  100  m  for  range,  r  <  47  km  and  r  >  52  km.  Ax  =  10  m  for  47 
km  <  r  <  50  km;  and  Ax  =  1  m  for  50  km  <  r  <  52  km  is  available  in  [Ref  1]  and  the 
Millington  approximation  in  [Ref  13].  We  computed  our  numerical  solution  based  on  the 
geometry  described  in  [Ref.  1].  In  view  of  the  large  distance,  the  asymptotic  form  of  Hankel 
function  was  used  in  Eq.  (76)  to  compute  the  magnetic  field  at  the  initial  range 
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where  R  =  |  kjix-x ^  \ .  The  path  loss  is  computed  as 
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and  plotted  in  Fig.  36.  We  see  that  there  is  excellent  agreement  in  the  path  loss  (attenuation 
factor)  between  our  computed  results  and  the  those  computed  using  the  IFDG  method  [Ref 
1].  We  also  noticed  that  the  general  behavior  of  our  numerical  solution  is  identical  to  the 
prediction  given  in  [Ref  1]. 


K.  PROPAGATION  OVER  A  LOSSY  VALLEY 

We  also  investigate  the  capability  of  the  parabolic  equation  formulation  to  predict 
radiowave  propagation  over  a  lossy  valley.  Consider  a  vertical  source  radiating  on  a  smooth 
horizontal  ground  at  a  distance  30  km  fi-om  the  edge  of  a  100-m-deep  valley  that  is  symmetric 
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about  the  midpoint  of  its  horizontal  bottom,  which  is  1  km  long.  The  slope  of  each  wall  is 
±1/10,  so  that  the  total  horizontal  length  of  the  terrain  irregularity  is  3  km.  The  ground  is 
characterized  by  =  4,  a  =  0.001  mho/m.  Results  using  the  IFDG  prediction  and  the 
WAGNER  method  [Ref  12]  for  receiver  on  the  ground  are  given  in  [Ref  1].  In  [Ref  1],  the 
IFDG  results  were  obtained  using  x,„,  =  25  km,  =  Im,  Ax  =  20  m  for  all  range,  r,  except 
for  30  km  <r  <35  km  for  which  Ax  =  1  m.  Again,  our  computational  geometry  is  similar  to 
the  IFDG  geometry.  The  initial  field  is  computed  using  Eq.  (76)  with  the  asymptotic  form  for 
Hankel  function  given  in  Eq.  (79).  The  path  loss  for  propagation  at  a  frequency  of  10  MHz 
is  computed  using  Eq.  (80)  and  plotted  in  Fig.  37.  Although  there  is  a  slight  deviation  in  the 
path  loss  between  our  numerical  results  and  the  IFDG  method  [Ref  1],  we  notice  that  the 
general  trend  of  the  path  loss  between  our  method  and  the  IFDG  method  remains  identical. 
The  deviation  in  the  magnitude  of  the  path  loss  could  be  due  to  the  fact  that  [Ref  12]  has 
considered  an  azimuthally  symmetrical  profile  while  our  profile  is  2-dimensional  and  laterally 
symmetrical. 
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Figure  36.  Path  loss  over  a  land-sea-land  path,/=  59.7  MHz,  transmitter  at  5m,  receiver 
at  5m,  vertical  polarization. 
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Figure  37.  Path  loss  for  propagation  over  a  valley,/=  10  MHz,  transmitter  and  receiver 
on  the  ground. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  thesis,  a  numerically  efficient  method  to  model  tropospheric  radiowave 
propagation  over  irregular  terrain  was  implemented  and  tested.  The  parabolic  equation 
method  was  used  for  predicting  radiowave  propagation.  A  tropospheric  boundary  condition 
was  used  to  truncate  the  computational  domain  at  the  top,  while  an  impedance  boundary 
condition  was  used  on  the  uneven  terrain  to  characterize  the  ground.  Because  the  parabolic 
equation  method  is  a  full-wave  method,  it  includes  all  aspects  of  wave  propagation  such  as 
forward  reflection,  reflection,  diflfiection,  and  surface  wave  propagation.  However,  it  ignores 
back-scattering.  Since  the  parabolic  equation  method  models  wave  propagation  only  in  the 
forward  direction,  it  allows  for  a  rapid  solution  of  the  fields  by  way  of  marching  along  in 
range  starting  from  an  initial  range.  A  great  advantage  of  the  parabolic  equation  method 
compared  to  the  commonly  used  ray  method  is  that  it  is  valid  in  the  shadow  region  where  the 
latter  method  completely  breaks  down.  Furthermore,  it  appears  to  be  the  only  practical 
method  for  predicting  propagation  over  long  ranges  (thousands  of  wavelengths)  and  over  a 
wideband  (HF  through  SHF). 

The  main  advancement  made  in  this  thesis  is  the  use  of  a  non-rectangular  mesh  in  a 
curvilinear  coordinate  system  to  model  the  uneven  terrain.  A  body  fitted  coordinate  is 
generated  based  on  the  specification  of  the  ground  profile,  and  the  parabolic  PDE,  together 
with  the  boundary  conditions,  are  cast  in  a  curvilinear  coordinate  system.  This  not  only  makes 
the  method  more  efficient,  but  also  permits  more  accurate  imposition  and  modeling  of  the 
boundary  conditions.  We  used  the  Crank-Nicolson  implicit  scheme  in  the  computational 
domain  to  solve  the  parabolic  PDE.  For  a  cartesian  grid,  this  method  is  convergent  and  stable 
for  all  finite  values  of  step  size.  In  this  method,  the  parabolic  PDE  is  considered  as  being 
satisfied  at  the  mid-point  of  the  computational  grid.  Central  difference  approximations,  which 
are  second-order  accurate,  are  used  in  the  interior  as  well  as  at  the  boundaries. 

Numerical  results  to  predict  radiowave  propagation  over  various  obstacles  were 
computed  and  validated  with  the  results  available  in  the  open  literature.  These  include  both 
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a  PEC  plane  and  impedance  ground  with  differing  permittivity  and  conductivity.  Excellent 
agreement  was  observed  and  demonstrated  for  the  following  cases  :  a  PEC  circular  boss, 
lossy  Gaussian  hill,  PEC  triangular  hill,  loaded  and  unloaded  cliff  edges,  mixed  path,  and  a 
lossy  valley.  Three-dimensional  obstacles  are  outside  the  scope  of  this  thesis. 

The  model  presented  in  this  thesis  can  be  applied  to  problems  in  short  range 
communications  or  radar  target  detection.  At  HFATIF  frequencies,  the  antennas  are 
electrically  close  to  the  ground  leading  to  significant  perturbation  of  the  received  signals. 
Assessment  of  propagation  of  HFA^HF  signals  over  inhomogeneous  and  irregular  terrain  is 
highly  desirable.  In  the  case  of  a  radar  system,  signals  received  from  a  target  depend  on  the 
direct  ray  and  a  ground  reflected  ray.  The  path  taken  by  the  ground  reflected  ray  depends  on 
the  terrain  roughness  and  ground  constants.  Since  target  detection  is  based  on  the  composite 
signal,  it  is  important  to  assess  ground  conditions  and  atmospheric  factors.  The  numerical 
model  developed  in  this  thesis  provides  a  useful  tool  for  predicting  the  performance  of  radar 
or  communication  links  before  commencing  to  design,  develop,  and  deploy  the  system  in  the 
field.  From  the  model  developed,  one  can  predict  propagation  factors  such  as  path  loss, 
coverage  diagrams,  and  the  perturbation  of  antenna  radiation  patterns  over  realistic  terrains 
and  ranges. 

This  thesis  represents  the  beginning  of  an  effort  to  predict  radiowave  propagation 
using  the  parabolic  equation  method.  In  its  standard  form,  the  accuracy  of  the  method  is 
limited  to  waves  traveling  within  ±15°  from  the  horizontal.  It  is  recommended  that  follow-on 
work  should  use  a  wide-angle  parabolic  equation  method  to  accommodate  waves  traveling 
over  wider  angles  and  in  an  inhomogeneous  atmosphere. 
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APPENDIX  A.  MATLAB  SOURCE  LISTINGS 


The  formulation  described  in  Chapter  II  and  the  numerical  procedure  described  in 
Chapter  III  were  implemented  using  MATLAB  sofware.  The  routines  that  generated  the 
results  in  Chapter  IV  are  listed  below. 


A.  EXECUTIVE  ROUTINE 

%  Filename :  main.m 

%  Title  :  Main  Program  or  Executive  Program 

%  Revision  No.  :  Rl.O 

%  Dale  of  Last  Revision  : 

10  Mar  95 

%  Comments : 

%  This  is  the  main  program  for  computing  tlie  Parabolic 

%  Equation  Solution  of  Radiowave  I’ropagation  in  an 

%  Inhomogenoiis  Atmosphere  and  (Jvcr  irregular  Teirain. 

% 

%  Input  Variables  :  none 

% 

%  Output  Variables  : 

%  xphy 

:  x-coordinate  of  physical  domain 

%  yphy 

:  y-coordinate  of  physical  domain 

%  U 

:  computed  Hz-Field 

%  wavelength 

:  wavelength  (frequency/3e8) 

%  time 

:  total  CPU  time  used 

% 

%  Associated  Matlab  Files  : 

%  Called  by 

:  none 

%  Subroutines 

:  data 

% 

dinp 

% 

cmet 

% 

hnd 

% 

epde 

% 

%  Associated  Functions 

:  none 

% 

clear  all; 

t=cputime; 

%  start  CPU  clock 

data 

%  input  data  file 

dinp 

%  compute  the  physical  domain 

cmet 

%  compute  the  metric 

hfld 

%  compute  Mz-Field  on  initial  data  line 

epde 

%  march  and  compute  Hz-Field 

time=cputime-t; 

save  output  xphy  yphy  N 

M  U  time  wavelengtlv. 

B.  PHYSICAL  DOMAIN  GRID 


%  Filename  :  dinp.m 

%  Title  :  Computation  of  Physical  Domain 
%  Revision  No.  :  Rl.O 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

%  This  subroutine  computes  the  entire  computational 

%  domain  based  on  the  user  input  data  for  the  physical 

%  ground  plane. 


:  x-coordinate  of  physical  domain 
:  y-coordinate  of  physical  domain 


%  Input  Variables : 

^0  >™  :  x-coordinate  of  ground  plane 

:  y-coordinate  of  ground  plane 

% 

%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

% 

%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  none 
% 

nseg=length(xin)- 1 ;  %  total  number  of  segment 

% 

%%%Vo%  Compute  horizontal  axis  of  the  physical  domain  %%%%% 
% 


%  total  number  of  segment 


st=I; 

for  i=  1  :nseg; 

npts(i)=round(sqrt((xin(i+ 1  )-xin(i))^2+(yin(i+ 1  )-yin(i))"^2)/xdel); 

xphy(st)=^(i); 

yphy(st)=yin(i); 

for  j=l  :npts(i)-l; 

xphy  (st+j  )=xin(i)+j  *(xin(i+ 1  )-xin(i))/npts(i); 
yphy(  1  ,st-Kj)==7in(i)+j*(yin(i+l  )-yin(i))/np^^ 
end; 

st=st+npts(i); 

end; 

xphy(st)=xin(length(xin));  %  include  the  last  point  for  x 
yphy(st)=yin(length(yin));  %  include  the  last  point  for  y 
% 


%%%%%  Compute  vertical  axis  of  the  physical  domain  %%%%% 
% 


M=length(xphy)- 1 ; 
N=round(yup/ydel); 
for  i=l:N; 


yphy(i-i-l  ,:)=(  1  -i/N).*yphy(  1  ,:)+(i/N)*yup; 
end; 
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C.  CURVILINEAR  COORDINATE  SYSTEM  METRIC 


%  Filename :  cmet.m 

%  Title  :  Generation  of  Curvilinear  Coordinate  System  Metric 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

%  This  subroutine  computes  the  various  metrics  or  derivatives 
%  (dx_xi,  dx_eta,  dy_xi,  dy_eta,  ddy_eta)  at  any  interior 

%  points.  In  the  Crank  Nicolson  implicit  scheme,  the  metrics 
%  are  needed  at  the  midpoint  with  respect  to  xi,  i.e.,  at 
%  xi=xi+del_xi/2. 

% 


%  Input  Variables : 
%  xphy 

%  yphy 

% 


:  x-coordinate  of  physical  domain 
:  y-coordinate  of  physical  domain 


%  Output  Variables : 
%  dx_xi 

%  dx_eta 

%  dy__xi 

%  <fy_eta 

%  ddy_eta 

% 


:  1  St  order  x-derivative  w.r.t.  xi 
:  1st  order  x-derivative  w.r.t.  eta 
:  1st  order  y-derivative  w.r.t.  xi 
:  1  st  order  y-derivative  w.r.t.  eta 
:  2nd  order  y-derivative  w.r.t.  eta 


%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  none 
% 

msize=ones(size(yphy(:,l ))); 
dx_xi=msize*(xphy(2:M+l)-xphy(l  :M)); 
dx_eta=0; 
for  i=l:N+l; 

dy_xi(i,l  :M)=(l-(i-l)/N)*(yphy(i>2:M+l)-yphy(i,l  :M)); 
end; 

dy_eta=msize*((yphy(N+l ,  1  )-(yphy(  1 ,2:M+1  )+yphy(  1 , 1  :M))./2)./N); 
ddy_eta=0; 


D.  MAGNETIC  FIELD  ON  INITIAL  DATA  LINE 

%  Filename :  hfld.m 

%  Title  :  Generation  of  Hz-Field  on  Initial  Data  Line 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 

%  This  subroutine  generates  the  magneic  field  (Hz-Field) 

%  at  the  initial  data  line  for  a  line  source  over 
%  an  impedance  plane. 

%  (Note  :  Never  place  the  initial  data  line  at  the 
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%  same  location  as  the  transmitter,  else  the  field 
%  generated  would  be  zero). 

% 

%  Input  Variables : 

%  yphy  :  y-coordinate  of  physical  domain 

% 

%  Output  Variables : 

^  Hz  :  Magnetic  Field  @  initial  data  line 

% 

%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  hankel2 
% 

yinit=[yphy(l  :N,1)'  yphy(N+l,l):ydel:ytop]; 

Ro=sqrt((xphy(  1  )-x0).''2+(yinit-y0).''2)*ko; 

% 

%%%%%  Generate  Hz-Field  using  Fredholm  Equation  %%%%% 

% 

hr2=hankel2(0  Jlo);  %  load  Hankel  function  of  2nd  order 

H=2  *  (( 1  j  *ko)/4) .  *((xphy(  1  )-xO)./Ro).  *hr2 ; 

Hz=H.'; 

% 

%  Note:  is  for  Nonconjugated  Transpose.  Using  only 

%  will  result  in  conjugated  transpose  (i.e.  the  sign 
%  of  the  complex  term  is  inverted) 


E.  FORMULATION  OF  PARABOLIC  EQUATION  MATRIX 

%  Filename  :  epde.m 

%  Title  :  Matrix  Equation  for  the  Parabolic  System  Equations 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 

%  This  program  defines  the  PDE  for  the  matrix  equations 
%  and  performs  the  numerical  implementation  using  the 
%  Crank-Nicolson  Implicit  Scheme  (Discretization). 

%  The  applicability  of  the  PDE  equation  is  also 

%  extended  to  accommodate  the  derivative  boundaiy 
%  conditions  on  the  lower  and  upper  boundaries. 

% 

%  Input  Variables : 


% 

dx_xi 

:  1  St  order  x-derivative  w.r.t.  xi 

% 

dx_eta 

:  1st  order  x-derivative  w.r.t.  eta 

% 

dy_xi 

:  1  St  order  y-derivative  w.r.t.  xi 

% 

dy_eta 

:  1  St  order  y-derivative  w.r.t.  eta 

% 

ddy  eta 

:  2nd  order  y-derivative  w.r.t.  eta 

% 

Hz 

:  Magnetic  Field  @  initial  data  line 

% 
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:  computed  Hz-Field 


%  Output  Variables : 

%  U 
% 

%  Associated  Matlab  Files  : 


% 

Called  by 

;  main 

% 

Subroutines 

:  elbc 

% 

eubc 

% 

%  Associated  Functions 

:  none 

% 

%%%%%  Assume  constant  atmospheric  refractive  index  %%%%% 
n=l ; dn_x=0; dn_y=0; 
al=-(2/n)*dn_x;  a2=-(2/n)*dn_y; 
a  1  _star=(n''2- 1  -2j*(a  1 /ko))*(ko''2); 

% 

%%%%%  Set  up  the  matrix  coefiBcients,  (M  x  N)  for  b  1 ,  b2,  b3  %%%%% 
% 

templ=2j*ko-al; 
b  1  =dx_xi.  *a  1  _star./temp  1 ; 

b2=dx_xi./dy_eta.*((a2-ddy_eta./dy_eta.''2)./templ+dy_xi./dx_xi); 

b3=dx_xi./(temp  1  .*dy_eta.''2); 

alpha=0.5.*(b3-K).5*b2); 

beta=(b3-0.5*bl); 

ganmia=0.5.*(b3-0.5.*b2); 

% 

%%%%%  Initialization  of  Magnetic  Field  %%%%% 

% 

U(l;N+l,l)=Hz(l:N+l); 

% 

%%%%%  Compute  coefficient  of  Lower  BC  %%%%% 

% 

elbc 

% 

%%%%%  Compute  invariant  coefficient  of  Upper  BC  %%%%% 

% 

r(l  :M)=sqrt((8j*ko)/pi)./sqrt((xphy(2:M+l)-xphy(l  :M))./2); 

% 

%%%%%  Loading  the  matrix  %%%%% 

% 

for  p=l:M; 

A=spdiags([gamma(:,p)  -( 1  +beta(:,p))  alpha(:,p)],  -1:1,  N+1 ,  N+1 ); 
B=spdiags([-gamma(;,p)  -(l-beta(:,p))  -alpha(:,p)],  -1:1,  N+1,  N+1); 

A(  1 ,2)=alpha_p  ( 1  ,p) ;  A(  1 , 1  )=-(  1  +beta_p(  1  ,p)); 

B(  1 ,2)=-alpha_p(  I  ,p);  B(  1 , 1  )=-(  1  -beta_pp(  1  ,p)); 

% 

%%%%%  Compute  coefficient  of  Upper  BC  %%%%% 

% 

eubc; 

lambda=beta(N+l  ,p)+dy_eta(p)*2*alpha(N+l  ,p)*(r(p)); 
ganima_p=gamma(N+ 1  ,p)+alpha(N+ 1  ,p); 

A(N+ 1  ,N+ 1  )=  1  +lambda; 

A(N+1  ,N)=-gammaj); 
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B(N+ 1  ,N+ 1  )=(  1  -lambda); 

B(N+1  ,N)=gamma_p; 

B(N+1,N+2)=1; 

U(N+2,p)=4*alpha(N+ 1  ,p)*dy_eta(p)*(s(p)); 

U(1  :N+l.p+l)=(A\B)*U(l  :N+2,p); 

end; 


F.  LOWER  BOUNDARY  CONDITION 


%  Filename :  elbc.m 
%  Title  :  Lower  Boundary  Conditions 
%  Revision  No. :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 


% 

% 

% 

% 

% 

% 

% 

% 

% 


This  program  formulates  the  equations  for  the 
impedance  boundaiy  conditions  imposed  on  the  irregular 
lower  boundary.  Impedance  boundary  conditions  relates 
the  tangential  components  of  electric  and  magnetic 
fields  at  the  interface  of  two  media. 

The  surface  impedance  is  determined  from  the 
intrinsic  impedance  of  the  medium  by  considering 
plane  wave  reflections  from  the  interface. 


%  Input  Variables  : 
%  dx_xi 

%  dy_xi 

%  dy_eta 

%  alpha 

%  gamma 

% 

%  Output  Variables  : 
%  alphajp 

%  beta_j) 

%  betajp 

% 


;  1st  order  x-derivative  w.r.t.  xi 
:  1  st  order  y-derivative  w.r.t.  xi 
:  1  st  order  y-derivative  w.r.t.  eta 
:  coefficient  for  U(p,q+1)  term 
;  coefficient  for  U(p,q- 1 )  term 


:  coefficient  for  U(p,  1 )  term 
:  coefficient  ”(l+beta_P)”  for  U(p,0)  term 
:  coefficient "( 1  -bela_pp)"  for  U(p- 1 ,0)  term 


%  Associated  Matlab  Files  : 


%  Called  by  :  epde 

%  Subroutines  :  none 

% 


%  Associated  Functions  :  none 
% 


temp2=sqrt(dx_xi(  1 , 1  :M).^2+dy_xi(  1 , 1  :M).^2); 
ct=dx_xi(  1 , 1  :M)./temp2;  %  cos(theta) 

st=dy_xi(  1 , 1  :M)  ./temp2 ;  %  sin(theta) 

cl=-lj*ko*(n^2.*delta_v  -  st); 
alpha j)=alpha(  1  ,:)+gamma(  1 
temp3=2.*gamma(l  ,;).*dy__eta(l  ,:).*ct; 
beta_j)=beta(  1  ,:)-temp3 .  *(c  1  -2*st./dx_xi(  1 , 1  :M)); 
beta  j)p=beta(  1  ,:)-temp3 .  *(c  1  +2*st./dx_xi(  1 , 1  :M)); 
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G. 


UPPER  BOUNDARY  CONDITION 


%  Filename :  eubc.m 
%  Title  :  Upper  Boundary  Conditions 
%  Revision  No.  :  Rl.O 
%  Date  of  Last  Revision  :  1 0  Mar  95 
%  Comments : 

%  This  program  formulates  the  equations  for  the  upper 
%  boundary  conditions  at  the  tropospheric  boundary. 

%  To  truncate  the  computational  domain  at  the  upper 

%  boundary,  a  point  high  enough  where  the  atmosphere 

%  is  homogeneous  is  considered. 

% 

%  Input  Variables : 

%  dx__xi 

%  dy_xi 

%  dy_eta 

%  alpha 

%  gamma 

% 

%  Output  Variables : 

%  s 

% 


:  1st  order  x-derivative  w.r.t.  xi 
:  1  St  order  y-derivative  w.r.t.  xi 
:  1  st  order  y-derivative  w.r.t.  eta 
:  coefficient  for  U(p,q+1 )  term 
:  coefficient  for  U(p,q- 1 )  term 


;  coefficient  for  U(p,  1 )  term 


%  Associated  Matlab  Files  : 

%  Called  by  :  epde 

%  Subroutines  ;  none 


% 

%  Associated  Functions  :  fresnel 


% 

tempA=0;  tempB=0;  tempC=0; 
x=(xphy(p+l)-i-xphy(p))/2;  %  x(p-l/2) 
xd=(xphy(p+ 1  )-xphy(p))/2;  %  x(p- 1  /2)-x(p- 1 ) 

for  i=N+2:round(ytop/ydel+l); 
j=i-N-l; 

temp4(j)=fresnel(sqrt(ko/(pi*(x-xin(  1 )))).  *(yinit(i)-yinit(N+ 1 ))); 
temp5(j)=fresnel(sqrt(ko/(pi*(x-xin(  1  )))).*  (yinit(i- 1  )-yinit(N+l ))); 
temp  A=tempA+(Hz(i,  1  )-Hz(i- 1 , 1  ))^del.  *(temp4(j)-temp5(j)); 
end; 

% 

%%%%%  Compute  the  U(N+1  ,p-l)  term  in  s(p)  %%%%% 

% 

tempA=sqrt(2j)*tempA; 

tempC=sqrt(8j  *ko/pi)*U  (N*^  1  ,p)/sqrt(xd); 

ifp=l; 

s(p)=tempA+tempC; 

else; 

temp6=sqrt(x-xphy(2:p))H-sqrt(x-xphy(l  :p- 1 )); 
tempB=sum((U(N+ 1 ,2:p)-U(N+l ,  1  :p- 1  ))./temp6); 
s(p)=tempA-sqrt(8j*ko/pi)*tempB+tempC; 
end; 
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H. 


FRESNEL  INTEGRAL  FUNCTION 


%  Filename :  fresnel.m 
%  Title  :  The  Fresnel  Integral  Function 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

% 

%  Input  Variables : 

% 

%  Output  Variables  : 

% 

%  Associated  Matlab  Files  : 

% 

%  Associated  Functions  : 

% 

function  F=fresnel(limits); 
interval=0.001; 
if  limits==0; 

F=0; 

elseif  limits  <0; 
limits—  1  *iimits; 
tau=0:interval:liniits; 
tau— l.*tau; 

integral=exp(- 1  j*pi/2*(tau.  *tau)); 
F=trapz(tau,integral); 
else 

tau=0  :interval;limits; 
integral=exp(- 1  j*pi/2*(tau.*tau)); 
F=trapz(tau,integral) ; 
end; 


L  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  PEC  PLANE 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


Filename :  pec.m 

Title  :  Input  Data  File 

Revision  No.  :  Rl.O 

Date  of  Last  Revision  :  10  Mar  95 

Comments : 

This  is  the  input  data  file  for  computing  the  Parabolic 
Equation  Solution  of  Radiowave  Propagation  in  an 
Inhomogenous  Atmosphere  and  Over  Irregular  Terrain. 
It  contains  the  frequency  of  operation,  the  x  and  y 
discretization  size,  x-  and  y-coordinates  of  the  ground 
plane  domain,  and  the  upper  boundary. 

It  also  contains  the  surface  impedance  (Note  :  surface 
impedance  =  0  for  PEC). 

Input  Variables  :  none 
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%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

%  wavelength  :  wavelength  (frequency/3  e8) 

% 

%  Associated  Matlab  Files  : 


%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  none 


% 


freq=300e6; 
wavelength=3e8/freq; 
ko=2*pi*freq/3e8; 
xdel=0. 1  ^wavelength; 
ydel=0. 1  *wavelength; 


%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


% 

%%%%%  Define  Ground  Plane  Coordinates  %%%%% 


% 

xin(  1  )=5*  wavelength;  yin(  1  )=0; 

xin(2)=15*wavelength;  yin(2)=0; 

yup=20*wavelength;  %  y-upper  boundary 

ytop=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 

x0=0;  y0=0;  %  position  of  transmitter 

% 

%%%%®/o  Define  Surface  Impedance  %%%%% 

% 

delta_v=0;  %  surface  impedance  (PEC) 


J.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  PEC  KNIFE-EDGE 

%  Filename ;  kpec.m 

%  Title  :  Input  Data  File  for  PEC  Knife  Edge 
%  Revision  No. :  R 1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

%  wavelength  :  wavelength  (frequency/3e8) 

% 

%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  none 
% 
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freq=300e6; 
wavelength=3e8/freq; 
ko=2*pi*freq/3e8; 
xdel=0. 1  ^wavelength; 
ydel=0. 1  ^wavelength; 
% 


%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


%%%%%  Define  Ground  Plane  Coordinates  %%%%% 

% 

xin(  1  )=5*wavelength;  ym(  1  )=0; 

xin(2)=9.9*wavelength;  ym(2)=0; 

xin(3)=  1 0*wavelength;  yin(3)=2*wavelength; 

xin(4)=  1 0. 1  *wavelength;  yin(4)=2*wavelength; 

xin(5)=  1 0.2*wavelength;yin(5)=0; 

xin(6)=30*wavelength;  yin(6)=0; 

yup=22*wavelength;  %  y-upper  boundary 

ytop=5*wavelength+yup;  %  troposphere  boundaiy  (@  infinity) 

x0=0;  y0=0;  %  position  of  transmitter 

% 


%%%%%  Define  Surface  Impedance  %%%%% 

% 

delta_v=0;  %  surface  impedance  (PEC) 


K.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  CIRCULAR  BOSS  IN  A 
PEC  PLANE 

%  Filename :  boss.m 

%  Title  :  Input  Data  File  for  Semi-Circular  Boss  over  PEC  Plane 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

% 

%  Input  Variables  ;  none 
% 

%  Output  Variables  : 


% 

xphy 

:  x-coordinate  of  physical  domain 

% 

yphy 

:  y-coordinate  of  physical  domain 

% 

% 

wavelength 

:  wavelength  (frequency/3e8) 

%  Associated  Matlab  Files  : 
%  Called  by 

%  Subroutines 

% 

%  Associated  Functions 
% 

freq=3e6; 

wavelength=3e8/jfreq; 
ko=2  *pi  *freq/3  e8 ; 
xdel=0.0 1  ^wavelength; 
ydel=0.05*wavelength; 

% 


:  main 
:  none 

:  none 

%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 
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%%%%%  Define  Ground  Plane  Coordinates  %%%%% 

% 

b=0 . 5  *  wavelength; 

A=wavelength;  B=0.8*wavelength; 
x_niid=-0.5*wavelength;0.0 1  *wavelength:0.5*wavelength; 
y_mid=sqrt(b^2-x_mid.^2); 

x_fet=-0.6*wavelength:0.0 1  *wavelength:-0.5 1  *wavelength; 
y_fiit=zeros(size(x_fiit)); 

x_bk=0.5 1  *wavelength:0.01  ^wavelength:  wavelength; 

y_bk=zeros(size(x_bk)); 

xin=[x_fiit  x_mid  x_bk]; 

%xin=xin-K).75*wavelength; 

yin=[y_fiit  y_niid  y_bk]; 

clear  x_fot;  clear  x_mid;  clear  x_bk; 

clear  y_fot;  clear  y_mid;  clear  y_bk; 

yup=20*wavelength;  %  y-upper  boundaiy 

ytop=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 

x0=-0.75*wavelength;  y0=0. 1  ^wavelength;  %  position  of  transmitter 

% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

delta_v=0;  %  surface  impedance  (PEC) 


L.  WPUT  DATA  FILE  FOR  PROPAGATION  OVER  LOSSY  IMPEDANCE 
PLANE 


%  Filename  :  imped.m 

%  Title  :  Input  Data  File  for  Lossy  Impedance  Plane 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 


% 


%  Input  Variables  :  none 
% 

%  Output  Variables : 


% 

xphy 

:  x-coordinate  of  physical  domain 

% 

yphy 

:  y-coordinate  of  physical  domain 

% 

wavelength 

:  wavelength  (frequency/3  e8) 

% 

%  Associated  Matlab  Files  ; 

% 

Called  by 

:  main 

% 

Subroutines 

:  none 

% 

%  Associated  Functions 

:  none 

% 

freq=le6; 

wavelength=3e8/freq; 
ko=2*pi*freq/3e8; 
xdel==0. 1  *wavelength; 
ydel=0. 1  *  wavelength; 


%  frequency  of  operation 
.  %  wavelength 

%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 
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% 

%%%%%  Define  Ground  Plane  Coordinates  %%%%% 

% 

xm(  1  )=0.5*wavelength;  yin(  1  )=0; 

xin(2)=20*wavelength;  yin(2)=0; 

yup=20*wavelength;  %  y-upper  boundary 

3^op=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 

x0=0;  y0=wavelength/300;  %  position  of  transmitter 

% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

er=  1 0;  %  relative  permittivity 

sigma=  1  Oe-3 ;  %  conductivity 

sigma_r=60*wavelength*sigma;  %  relative  conductivity 

delta_v=sqrt(l/(er-lj*sigma_r));%  surface  impedance  (vertical  pol) 


M.  INFVT  DATA  FILE  FOR  PROPAGATION  OVER  LOSSY  GAUSSIAN  HILL 


%  Filename :  gauss.m 

%  Title  :  Input  Data  File  for  Lossy  Gaussian  Hill 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

%  wavelength  :  wavelength  (frequency/3  e8) 

% 


%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 


%  Associated  Functions 
% 

freq=le6; 

wavelength=3e8/freq; 
ko=2*pi’*‘freq/3e8; 
xdel=0. 1/3*  wavelength; 
ydel=0. 1  *wavelength; 

% 


:  none 


%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


%%%%%  Define  Ground  Plane  Coordinates  %%%%% 


% 

Ht=  1 000;  c=5000;  w=3000;  %  constants  for  gaussian  hill 

xin=2000:xdel:24*wavelength;%  x-coordinate 
yin=Ht.*exp(-9.*(((xin-c)./w).^2));  %  y-coordinate 
yup=30*wavelength;  %  y-upper  boundary 

ytop=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 
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x0=0;  y0=wavelength/100;  %  position  of  transmitter 
% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

eT=  1 0;  %  relative  permittivity 

sigma=  1  Oe-3 ;  %  conductivity 

sigma_r=60*wavelength*sigma;  %  relative  conductivity 

delta_v=sqrt(l/(er-lj*sigma_r));%  surface  impedance  (vertical  pol) 


N.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  PEC  ISOSCELES  HILL 


%  Filename :  isos.m 

%  Title  :  Input  Data  File  for  PEC  Isosceles  Hill 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

%  wavelength  :  wavelength  (frequency/3  e8) 

% 


%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 


% 

%  Associated  Functions 
% 

freq=9.6e9; 
wavelength=3  eS/freq; 
ko=2*pi*freq/3e8; 
xdel=0.08*wavelength; 
ydel=0.08*wavelength; 
% 


:  none 


%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


%%%%%  Define  Ground  Plane  Coordinates  %%%%% 


% 

xin(  I  )=5. 1 2*wavelength;  yin(  1  )=0; 

xin(2)=  10.24*  wavelength;  yin(2)=0; 

xin(3)=20.32*wavelength;  yin(3)=l  .232*wavelength; 

xin(4)=30.4*wavelength;  yin(4)=0; 

xin(5)=64*wavelength;  yin(5)=0; 

yup=20*wavelength;  %  y-upper  boundary 

ytop=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 

x0=0;  yOK);  %  position  of  transmitter 

% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

delta_v=0;  %  surface  impedance  (PEC) 
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O.  mPUT  DATA  FILE  FOR  PROPAGATION  OVER  PEC  CLIFF 


%  Filename :  cliff.m 
%  Title  :  Input  Data  File  for  PEC  Cliff 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 


% 

xphy 

:  x-coordinate  of  physical  domain 

% 

yphy 

:  y-coordinate  of  physical  domain 

% 

% 

wavelength 

:  wavelength  (frequency/3  e8) 

%  Associated  Matlab  Files  : 

%  Called  by 

%  Subroutines 

% 

%  Associated  Functions 
% 

freq=9.6e9; 
wavelength=3e8/freq; 
ko=2*pi*freq/3e8; 
xdel=0.08*wavelength; 
ydel=0. 08*  wavelength; 

% 

%%%%%  Define  Ground  Plane  Coordinates  %%%%% 

% 

xin(  1  )=5. 1 2*wavelength;  yin(  1  )=  1 .0464*  wavelength; 

xin(2)=  1 0.24*wavelength;  yin(2)=  1 ,0464*wavelength; 

xin(3)=10.34*wavelength;  yin(3)=0; 

xin(4)=80*wavelength;  yin(4)=0; 

yup=2 1  *wavelength;  %  y-upper  boundaiy 

ytop=5*wavelength+yup;  %  troposphere  boundary  (@  infinity) 

x0=0;  y0=yin(  1 );  %  position  of  transmitter 

% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

delta_v=0;  %  surface  impedance  (PEC) 


;  main 
:  none 

:  none 

%  fi-equency  of  operation 
%  wavelength 
%  fi‘ee  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


P.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  DIELECTRIC  COATED 
CLIFF 

%  Filename ;  diele.m 

%  Title  :  Input  Data  File  for  Dielectric  coated  Cliff 
%  Revision  No.  :  Rl.O 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments  : 
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% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 

%  xphy  :  x-coordinate  of  physical  domain 

%  yphy  :  y-coordinate  of  physical  domain 

%  wavelength  :  wavelength  (frequency/3e8) 

% 


%  Associated  Matlab  Files  : 

%  Called  by  :  main 

%  Subroutines  :  none 

% 

%  Associated  Functions  :  none 


% 

freq=9.6e9; 
wavelength=3  eS/freq; 
ko=2*pi*freq/3e8; 
xdel=0.08*wavelength; 
ydel=0.08*  wavelength; 
% 


%  frequency  of  operation 
%  wavelength 
%  free  space  wave  number 
%  size  of  delta-x 
%  size  of  delta-y 


%%%%%  Define  Ground  Plane  Coordinates  %%%%% 
% 


xin(  1  )=0.5  ^wavelength;  yin(  1  )=  1 .0464*wavelength; 

xin(2)=  1 0.24*wavelength;  yin(2)=  1 .0464*wavelength; 

xin(3)=10.34*wavelength;  yin(3)=0; 

xm(4)=80*wavelength;  yin(4)=0; 

yup=2 1  *wavelength;  %  y-upper  boundary 

ytop=5*wavelength+yup;  %  troposphere  boundaiy  (@  infinity) 

x0=0;  yO=yin(  1 );  %  position  of  transmitter 

% 


%%%%%  Define  Surface  Impedance  %%%%% 

% 

delta_v=0. 1 7 1  j;  %  surface  impedance 


Q.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  MIXED  PATH 

%  Filename :  mix.m 

%  Title  :  Input  Data  File  for  Land-Sea-Land 
%  Revision  No.  :  R1 .0 
%  Date  of  Last  Revision  :  10  Mar  95 
%  Comments : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables  : 


% 

xphy 

:  x-coordinate  of  physical  domain 

% 

yphy 

:  y-coordinate  of  physical  domain 

% 

% 

wavelength 

:  wavelength  (frequency/3e8) 

%  Associated  Matlab  Files  : 
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%  Called  by 

%  Subroutines 

% 

%  Associated  Functions 
% 

freq=59.7e6; 

wavelength=3e8/freq; 

ko=2*pi*freq/3e8; 

ydel=0.2; 

% 


:  main 
;  none 

:  none 


%  wavelength 
%  Free  space  wave  number 
%  y-delta  size 


%%%%%  Define  Ground  Plane  Coordinates  %%%%% 
% 


xin(l)=40000;  yin(l)=0;  xdel(l)=100;  impz(l)=sqrt(l/(15-9.045j)); 
xin(2)=47000;  ym(2)=0;  xdel(2)=10;  impz(2)=sqrt(l/(80-1206j)); 
xm(3)=50000;  yin(3)=0;  xdel(3)=l ;  impz(3)=sqrt(l/(15-9.045j)); 
xin(4)=52000;  yin(4)=0;  xdel(4)=100;  impz(4)=sqrt(l/(15-9.045j)); 
xin(5)=60000;  yin(5)=0;  impz(5)=sqrt(l/(15-9.045j)); 
yup=10; 

ytop=5*wavelength4yup; 

x0=0;y0=5; 


R.  INPUT  DATA  FILE  FOR  PROPAGATION  OVER  LOSSY  VALLEY 

%  Filename :  valley.m 
%  Title  :  Input  Data  File  for  Valley 
%  Revision  No.  ;  R0.2 
%  Date  of  Last  Revision  :  20  Feb  95 
%  Comments  : 

% 

%  Input  Variables  :  none 
% 

%  Output  Variables : 


%  xphy 

:  x-coordinate  of  physical  domain 

%  yphy 

:  y-coordinate  of  physical  domain 

%  wavelength 

:  wavelength  (frequency/3  e8) 

% 

%  Associated  Matlab  Files 

%  Called  by 

:  main 

%  Subroutines 

:  none 

% 

%  Associated  Functions 

:  none 

% 

freq=10e6; 

%  frequency  of  operation 

wavelength=3e8/freq; 

%  wavelength 

ko=2  *pi  *freq/3  e8 ; 

%  free  space  wave  number 

ydel=0. 1  ^wavelength; 

%  size  of  delta-y 

% 

%%%%%  Define  Ground  Plane  Coordinates  %%%%% 
% 

xin(  I  )=28000;  yin(  1  )=  1 00;  xdel(  I  )=300; 
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xin(2)=30000;  ym(2)=100;  xdel(2)=3; 
xm(3)=31000;yin(3)=0;  xdel(3)=3; 
xiii(4)=32000;  yin(4)=0;  xdel(4)=3; 
xin(5)=33000;  yin(5)=100;  xdel(5)=3; 
xm(6)=35000;  ym(6)=100;  xdel(6)=30; ' 
xiii(7)=50000;  yin(7)=100; 
yup=50+yin(l); 
ytop=5  *wavelength+yup; 
xO=<);yO=yin(l); 

% 

%%%%%  Define  Surface  Impedance  %%%%% 

% 

ei=4;  %  relative  permittivity 

sigma=  1  e-3 ;  %  conductivity 

sigma_r=60*wavelength*sigma;  %  relative  conductivity 

delta_v=sqrt(l/(er-lj*sigma_r));%  surface  impedance  (vertical  pol) 
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